home *** CD-ROM | disk | FTP | other *** search
/ CD Ware Multimedia 1995 May / cd Ware (Juegos) Epimundo.iso / DOS / C / CMATH.ZIP / LOG.C < prev    next >
Encoding:
C/C++ Source or Header  |  1989-03-21  |  6.5 KB  |  323 lines
  1. /*                            log.c
  2.  *
  3.  *    Natural logarithm
  4.  *
  5.  *
  6.  *
  7.  * SYNOPSIS:
  8.  *
  9.  * double x, y, log();
  10.  *
  11.  * y = log( x );
  12.  *
  13.  *
  14.  *
  15.  * DESCRIPTION:
  16.  *
  17.  * Returns the base e (2.718...) logarithm of x.
  18.  *
  19.  * The argument is separated into its exponent and fractional
  20.  * parts.  If the exponent is between -1 and +1, the logarithm
  21.  * of the fraction is approximated by
  22.  *
  23.  *     log(1+x) = x - 0.5 x**2 + x**3 P(x)/Q(x).
  24.  *
  25.  * Otherwise, setting  z = 2(x-1)/x+1),
  26.  * 
  27.  *     log(x) = z + z**3 P(z)/Q(z).
  28.  *
  29.  *
  30.  *
  31.  * ACCURACY:
  32.  *
  33.  *                      Relative error:
  34.  * arithmetic   domain     # trials      peak         rms
  35.  *    IEEE      0.5, 2.0    35000       1.8e-16     5.6e-17
  36.  *    IEEE      1, MAXNUM   10000       1.8e-16     4.8e-17
  37.  *    DEC       0.5, 2.0    20000       2.0e-17     7.0e-18
  38.  *    DEC       1, MAXNUM   17700       1.9e-17     6.1e-18
  39.  *
  40.  * In the tests over the interval [1, MAXNUM], the logarithms
  41.  * of the random arguments were uniformly distributed over
  42.  * [0, MAXLOG].
  43.  *
  44.  * ERROR MESSAGES:
  45.  *
  46.  * log singularity:  x = 0; returns MINLOG
  47.  * log domain:       x < 0; returns MINLOG
  48.  */
  49.  
  50. /*
  51. Cephes Math Library Release 2.1:  December, 1988
  52. Copyright 1984, 1987, 1988 by Stephen L. Moshier
  53. Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  54. */
  55.  
  56. #include "mconf.h"
  57. static char fname[] = {"log"};
  58.  
  59. /* Coefficients for log(1+x) = x - x**2/2 + x**3 P(x)/Q(x)
  60.  * 1/sqrt(2) <= x < sqrt(2)
  61.  */
  62. #ifdef UNK
  63. static double P[] = {
  64.   4.58482948458143443514E-5,
  65.   4.98531067254050724270E-1,
  66.   6.56312093769992875930E0,
  67.   2.97877425097986925891E1,
  68.   6.06127134467767258030E1,
  69.   5.67349287391754285487E1,
  70.   1.98892446572874072159E1
  71. };
  72. static double Q[] = {
  73. /* 1.00000000000000000000E0, */
  74.   1.50314182634250003249E1,
  75.   8.27410449222435217021E1,
  76.   2.20664384982121929218E2,
  77.   3.07254189979530058263E2,
  78.   2.14955586696422947765E2,
  79.   5.96677339718622216300E1
  80. };
  81. #endif
  82.  
  83. #ifdef DEC
  84. static short P[] = {
  85. 0034500,0046473,0051374,0135174,
  86. 0037777,0037566,0145712,0150321,
  87. 0040722,0002426,0031543,0123107,
  88. 0041356,0046513,0170752,0004346,
  89. 0041562,0071553,0023536,0163343,
  90. 0041542,0170221,0024316,0114216,
  91. 0041237,0016454,0046611,0104602
  92. };
  93. static short Q[] = {
  94. /*0040200,0000000,0000000,0000000,*/
  95. 0041160,0100260,0067736,0102424,
  96. 0041645,0075552,0036563,0147072,
  97. 0042134,0125025,0021132,0025320,
  98. 0042231,0120211,0046030,0103271,
  99. 0042126,0172241,0052151,0120426,
  100. 0041556,0125702,0072116,0047103
  101. };
  102. #endif
  103.  
  104. #ifdef IBMPC
  105. static short P[] = {
  106. 0x974f,0x6a5f,0x09a7,0x3f08,
  107. 0x5a1a,0xd979,0xe7ee,0x3fdf,
  108. 0x74c9,0xc66c,0x40a2,0x401a,
  109. 0x411d,0x7e3d,0xc9a9,0x403d,
  110. 0xdcdc,0x64eb,0x4e6d,0x404e,
  111. 0xd312,0x2519,0x5e12,0x404c,
  112. 0x3130,0x89b1,0xe3a5,0x4033
  113. };
  114. static short Q[] = {
  115. /*0x0000,0x0000,0x0000,0x3ff0,*/
  116. 0xd0a2,0x0dfb,0x1016,0x402e,
  117. 0x79c7,0x47ae,0xaf6d,0x4054,
  118. 0x455a,0xa44b,0x9542,0x406b,
  119. 0x10d7,0x2983,0x3411,0x4073,
  120. 0x3423,0x2a8d,0xde94,0x406a,
  121. 0xc9c8,0x4e89,0xd578,0x404d
  122. };
  123. #endif
  124.  
  125. #ifdef MIEEE
  126. static short P[] = {
  127. 0x3f08,0x09a7,0x6a5f,0x974f,
  128. 0x3fdf,0xe7ee,0xd979,0x5a1a,
  129. 0x401a,0x40a2,0xc66c,0x74c9,
  130. 0x403d,0xc9a9,0x7e3d,0x411d,
  131. 0x404e,0x4e6d,0x64eb,0xdcdc,
  132. 0x404c,0x5e12,0x2519,0xd312,
  133. 0x4033,0xe3a5,0x89b1,0x3130
  134. };
  135. static short Q[] = {
  136. 0x402e,0x1016,0x0dfb,0xd0a2,
  137. 0x4054,0xaf6d,0x47ae,0x79c7,
  138. 0x406b,0x9542,0xa44b,0x455a,
  139. 0x4073,0x3411,0x2983,0x10d7,
  140. 0x406a,0xde94,0x2a8d,0x3423,
  141. 0x404d,0xd578,0x4e89,0xc9c8
  142. };
  143. #endif
  144.  
  145. /* Coefficients for log(x) = z + z**3 P(z)/Q(z),
  146.  * where z = 2(x-1)/(x+1)
  147.  * 1/sqrt(2) <= x < sqrt(2)
  148.  */
  149.  
  150. #ifdef UNK
  151. static double R[3] = {
  152. -7.89580278884799154124E-1,
  153.  1.63866645699558079767E1,
  154. -6.41409952958715622951E1,
  155. };
  156. static double S[3] = {
  157. /* 1.00000000000000000000E0,*/
  158. -3.56722798256324312549E1,
  159.  3.12093766372244180303E2,
  160. -7.69691943550460008604E2,
  161. };
  162. #endif
  163. #ifdef DEC
  164. static short R[12] = {
  165. 0140112,0020756,0161540,0072035,
  166. 0041203,0013743,0114023,0155527,
  167. 0141600,0044060,0104421,0050400,
  168. };
  169. static short S[12] = {
  170. /*0040200,0000000,0000000,0000000,*/
  171. 0141416,0130152,0017543,0064122,
  172. 0042234,0006000,0104527,0020155,
  173. 0142500,0066110,0146631,0174731,
  174. };
  175. #endif
  176. #ifdef IBMPC
  177. static short R[12] = {
  178. 0x0e84,0xdc6c,0x443d,0xbfe9,
  179. 0x7b6b,0x7302,0x62fc,0x4030,
  180. 0x2a20,0x1122,0x0906,0xc050,
  181. };
  182. static short S[12] = {
  183. /*0x0000,0x0000,0x0000,0x3ff0,*/
  184. 0x6d0a,0x43ec,0xd60d,0xc041,
  185. 0xe40e,0x112a,0x8180,0x4073,
  186. 0x3f3b,0x19b3,0x0d89,0xc088,
  187. };
  188. #endif
  189. #ifdef MIEEE
  190. static short R[12] = {
  191. 0xbfe9,0x443d,0xdc6c,0x0e84,
  192. 0x4030,0x62fc,0x7302,0x7b6b,
  193. 0xc050,0x0906,0x1122,0x2a20,
  194. };
  195. static short S[12] = {
  196. /*0x3ff0,0x0000,0x0000,0x0000,*/
  197. 0xc041,0xd60d,0x43ec,0x6d0a,
  198. 0x4073,0x8180,0x112a,0xe40e,
  199. 0xc088,0x0d89,0x19b3,0x3f3b,
  200. };
  201. #endif
  202.  
  203.  
  204. #define SQRTH 0.70710678118654752440
  205. extern double MINLOG;
  206.  
  207. double log(x)
  208. double x;
  209. {
  210. int e;
  211. short *q;
  212. double y, z;
  213. double frexp(), ldexp(), polevl(), p1evl();
  214.  
  215. /* Test for domain */
  216. if( x <= 0.0 )
  217.     {
  218.     if( x == 0.0 )
  219.         mtherr( fname, SING );
  220.     else
  221.         mtherr( fname, DOMAIN );
  222.     return( MINLOG );
  223.     }
  224.  
  225. /* separate mantissa from exponent */
  226.  
  227. #ifdef DEC
  228. q = (short *)&x;
  229. e = *q;            /* short containing exponent */
  230. e = ((e >> 7) & 0377) - 0200;    /* the exponent */
  231. *q &= 0177;    /* strip exponent from x */
  232. *q |= 040000;    /* x now between 0.5 and 1 */
  233. #endif
  234.  
  235. /* Note, frexp is used so that denormal numbers
  236.  * will be handled properly.
  237.  */
  238. #ifdef IBMPC
  239. x = frexp( x, &e );
  240. /*
  241. q = (short *)&x;
  242. q += 3;
  243. e = *q;
  244. e = ((e >> 4) & 0x0fff) - 0x3fe;
  245. *q &= 0x0f;
  246. *q |= 0x3fe0;
  247. */
  248. #endif
  249.  
  250. /* Equivalent C language standard library function: */
  251. #ifdef UNK
  252. x = frexp( x, &e );
  253. #endif
  254.  
  255. #ifdef MIEEE
  256. x = frexp( x, &e );
  257. #endif
  258.  
  259.  
  260.  
  261. /* logarithm using log(x) = z + z**3 P(z)/Q(z),
  262.  * where z = 2(x-1)/x+1)
  263.  */
  264.  
  265. if( (e > 2) || (e < -2) )
  266. {
  267. if( x < SQRTH )
  268.     { /* 2( 2x-1 )/( 2x+1 ) */
  269.     e -= 1;
  270.     z = x - 0.5;
  271.     y = 0.5 * z + 0.5;
  272.     }    
  273. else
  274.     { /*  2 (x-1)/(x+1)   */
  275.     z = x - 0.5;
  276.     z -= 0.5;
  277.     y = 0.5 * x  + 0.5;
  278.     }
  279.  
  280. x = z / y;
  281.  
  282.  
  283. /* rational form */
  284. z = x*x;
  285. z = x + x * ( z * polevl( z, R, 2 ) / p1evl( z, S, 3 ) );
  286.  
  287. goto ldone;
  288. }
  289.  
  290.  
  291.  
  292. /* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */
  293.  
  294. if( x < SQRTH )
  295.     {
  296.     e -= 1;
  297.     x = ldexp( x, 1 ) - 1.0; /*  2x - 1  */
  298.     }    
  299. else
  300.     {
  301.     x = x - 1.0;
  302.     }
  303.  
  304.  
  305. /* rational form */
  306. z = x*x;
  307. y = x * ( z * polevl( x, P, 6 ) / p1evl( x, Q, 6 ) );
  308. y = y - ldexp( z, -1 );   /*  y - 0.5 * z  */
  309. z = x + y;
  310.  
  311. ldone:
  312.  
  313. /* recombine with exponent term */
  314. if( e != 0 )
  315.     {
  316.     y = e;
  317.     z = z - y * 2.121944400546905827679e-4;
  318.     z = z + y * 0.693359375;
  319.     }
  320.  
  321. return( z );
  322. }
  323.